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ABSTRACT 


Recent research by narrison on the stability of parallel flows re- 
sulted in a successful solution of plane Poiseuille flow but produced 
unexplained anomalies for pipe flow. The purpose of the research in 
this paper was to find and correct errors in Harrison's initial analysis 
of the pipe flow problem. 

A minor error in Harrison's numerical method was corrected. More- 
over, it proved necessary to reanalyze the conditions on the axis of 
symmetry of the pipe. These changes finally made it possible to obtain 
reasonable results. 

Owing to time limitations, the number of solutions obtained using 
the corrected program was sufficient only to confirm its general vali- 
G@ity. however, the results obtained are significant in that they dis- 
close instabilities which are known to exist but which have not been 


accounted for in previous theoretical investigations. 
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TABLE OF SYMBOLS 


All quantities below are in dimensionless form, 


The partial derivatives with respect to r in cylindrical 


coordinates. 


Base of natural logarithms. Exponentiation is also de- 


noted by exp( ). 


Unit vectors along the x,r, and @ axes in cylindrical 


coordinates. 


Components of the velocity vector potential defined in 


equation (2-1). 


+V¥ 1 » the imaginary unit. Also used as an index in the 


finite difference mesh in Section IV. 


The number of interior points in the finite difference 
mesh of Section IV. Also used as the imaginary part of 


B, dee D=1N. 


Reynolds number based on the mean velocity and pipe ra- 


dius. 


Time 


T,t,P,9 


x,r,9 


Lo lelg 


Sep 


shorthand notation for commonly occurring groups of 
symbols defined in equations (2-18), (2-19), (3-2), 


and (3-3). 


Velocity vector of the perturbation flow defined by 


equation (2-3). 


Complex vector potential of perturbation velocity defined 


by equations (2-1) and (2-2). 
Cylindrical coordinates. 


Ay + 1A, Complex wave number of the perturbation in the 


x direction. 


iB, Complex wave number of the perturbation in the © di- 


I 


rection 
Vp + By Complex frequency of the perturbation. 


The vorticity transport eauation (2-5) expressed in ab- 


breviated notation as defined by equation (2-6). 


The components of Fin cylindrical coordinates defined by 


equation (2-6). 


Vorticity vector of the perturbation flow defined in 


equation (2-l). 


J 


{} 


Linear vector operator (nabla). 


Vector cross-product operator. 


Brackets enclosing a matrix. 


Brackets enclosing a column vector. 
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I. INTRODUCTION 


The research reported in this thesis was undertaken in an attempt 
to isolate sources of error causing obviously incorrect ccmputer so- 
lutions for the three-dimensional linearized vorticity transport equa- 
tion for pipe Poiseuille flow. Solutions obtained using the theory and 
program developed in Ref. 1 had indicated decreasing flow stability 
With decreasing Reynolds number, a result which is clearly inconsistent 
with theory and experimental data. 

Initial analysis of that program confirmed that it was basically 
correct as presented in Ref. 1, but the nature of the solutions ob- 
tained indicated one or more errors in the formulation of the problem 
to be solved. Analysis of this problem formulation was conducted in 
three phases. The first phase consisted of ensuring accurate expres- 
Sion of the vorticity transport equation in terms of a complex velocity 
vector potential. The second and most significant phase of the analy- 
sis required extensive examination of the conditions required to satis- 
fy the vorticity transport equation at the singular point on the axis 
of symmetry of the pipe. The third phase entailed establishing a fi- 
nite difference approximation of this equation and its boundary con- 
ditions for adaptation of the problem to a standard form for computer 


Solution. 


11 


II. TEE VORTICITY TRANSPORT EQUATION 


The governing equations for laminar flow of an incompressible fluid 
with constant viscosity are the continuity equation and the Navier- 
Stokes equation. Taking the curl (vx) of the latter equation and intro- 
ducing a perturbation velocity (v) and vorticity (@) leads to the 
linearized perturbation vorticity transport equation as derived in 
mocendix A of Ref. 1. 

Expressed in terms of the complex velocity vector potential, W, 


which has the form 


W(x,r,0,t) = (6F(r) + 6 G(r) + &H(r))e% CA) 
where 

ho] Ax 4IEO 4+4¥i (2-2) 
and 

v= vxi (2-3) 

G) = vxv P (2-l) 


the vorticity transport equation becomes three simultaneous fourth-order 


differential equations of the form 


pir DF D°F 
LY) J pig b « [N3] pa} + (TH) + ine) na b+ 
Dy DH DH 
DF F 0 
(IM) +¥LN,9) DG P + (Cio5 *¥LN)) G+ = <0 (2-5) 
DH H 0 


as derived in Appendix E of Ref. 1. 


Equations (2-5) may be expressed in the abbreviated form 


Wa 


Ee 0 
ele i este (2-6) 
is © 


where M appears to be a set of three coupled equations in the components 
of W. As shown in Appendix B of Ref. 1, equations (2-6) actually ex- 
press only two independent conditions. Therefore, an appropriate linear 


combination of [, and [' allows expression of equations (2-5) as a set 


Q 


of two equations in three unknowns. The linear combination 


Be Ne 3 (2-7) 
ty 


does not, in general, uncouple Pf but does reduce the highest order deri- 
vative of the component G(r) in equations (2-5) to second order. 

In a manner similar to that just described for, Appendix C of 
Ref. 1 iliustrates the redundancy of the three components of We This 
aliows arbitrary selection of one of these components as being uniformly 
zero for all r. As will be seen later, the maximum benefits from the 
consequences of the linear combination expressed by equation (2-7) are 
obtained if 

aor) = O ; (2-8) 

Incorporation of equations (2-7) and (2-8) in equation (2-5) results 
in expression of the vorticity transport equation in the form 


Weg 


(pte (°c | ; 
on py ance N° + (OMo} - sO8o) ma : 


,' DG ' G 0 
(CMJ -s£8,3) a + ((Mo] “ation |} = {-{ (2-9) 


where 


(2-10) 


[m,J = 10 -2A (2-11) 
rRe 
O B | 


(1,1) - iy at + na ie - }) (2-12) 
2 Re 2 
Dane 1 


pa 

[a 

Gea 

i 

© © 

CO ae 
D 


ot 2B 
fe “Re 
{ 
[uj J- [ua 2({at* ae - ‘ ‘ (2-13) 
3 r\ Ref 2 
r~ Re 1e 
(5 - ) 4 ) 
rRe ‘a r 2 
r r Re 
Hl 
LM) 7 (2-ih) 


2 (e{* * Si) ‘( afe [a] >) 
2 Re 2 ( 2 | 2 
x if ats 1 
Pie al (* =) (afk * 5] - i 
2 2 2 Re 2 
Tate 6 ie 1G 
(",) = 10 -A (2-15) 
oO 0 


[N,} = O -A (2-16) 
ie 
es 
is 


[No] = | 22AB ‘(27 ') (2-17) 
2 2 
1g i & 
ip By 
2 
is 


and 


a= + 8 (228) 


T = 2A(1 - a . (2-19) 


t 
R 


(D 


This statement of the vorticity transport equation is identical to 
that published in Ref. 1. Note that equation (2-8) makes it unneces- 
sary to include F(r) and its derivatives in equations (2-9). The co- 
efficients of this function and its derivatives are included in the 
arrays represented by equations (2-15) through (2-22) of Ref. 1 and 
have been verified. The reader's attention is called to a misprint 
ocurring in Ref. 1 concerning the development to this point. Equation 


(2-22) of Ref. 1 should be corrected to read as follows. 


[uo) = [Bt  2aB ‘(‘ - inh (2-20) 
a c 2 
I ie 
0 t -B 
2 
We 
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JIL. BOUNDARY AND SPECIAL CONDITIONS 


A. BOUNDARY CONDITIONS AT THE WALL 

The boundary conditions at the wall, r=1, are derived from the fact 
that the axial, radial, and angular components of the flow velocity at 
the wall are identically zero. These components are obtained respec- 
tively from the three components of the velocity vector resulting from 
equation (2-3). 

As shown in part two of Appendix F in Ref 1 for the case F(r)=0, 


these components imply the following boundary conditions at the wall. 


G(1) = 0 (3-12) 
Hil) = 0 (3-1b) 
DH 1} = 0 : (3-1c) 


B. CONDITIONS ON THE AXIS OF SYMMETRY 

{It is important to note that the line of points constituting the .~ 
axis of symmetry is not a boundary but rather a line consisting of an 
infinite number of singular points in the flow. 

It should be noted in passing that all functions are periodic with 
respect to @. Consequently, as shown in part one of Appendix G in Ref. 
1, the wave number of the perturbation, B, must be a pure imaginary 
number such that 

B=ni, Tn ORM ngecr serene (3-2) 
Hence references to B in this paper shall be understood to indicate a 
pure imaginary number. 

The primary condition to be imposed is that the vorticity transport 


equation be satisfied at the point r=O. Inspection of equations (2-9) 
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through (2-19) shows that, in general, these are not satisfied for 


@arect substitution of r=0. To resolve this situation, consider the 


following development. 


The coefficients of the functions G(r) and H(r) and their deriva- 


tives with respect to r uv to fourth order in equations (2-9) through 


(2-19) may be expressed in the following form. 
[4] eat 
[1] 32 Ue 
() = Cx hve" eee We ae n° 
04) = (have ely a = 
[%o) = [Nove +7 pole" eal [Moale® 
[71 = C0] 
cn See = 
[¥od = [o-2d/** * Vo-0] 


where 
[%-0) > | ° 
0 
(¥, 4) = 16 aa 
oe 
Re 


-lAB A (3 = 2B°) 
Re Re 


[Moo] = 


(3-3) 
(3-li) 
(3-5) 
(3-6) 
(3-7) 
(3-8) 
(3-9) 
(3-10) 


O-1) 


Cote) 


(3-13) 


oe 7 0 af i a) (3-1)) 


Re 
=ae 0 
Re 
: _ 2 
[5,0] = }9 -2aA (3-15) 
O 0 
2 
M = |LAB A (2B - 3 (3-16) 
W-3] Re Re 
BoB (BS = 1) 
Re Re 


1 - =) (3-17) 


(3-18) 
A_ pe) (3-19) 
Re 
B (Be + 1) 


onc % - As) 6) (3-20) 
Re 


+ ao ABBA - 2 
Re Re 


(3-21) 


(Np,0) = 24! (3+22) 
is 0 


(%,_9) = a (3-23) 
0 

EyojJ= [oO A (3-2h) 
B 

[No_o| = es A(1 - B _ (3~25) 


ieee = | O °) (3-26) 
i «0 


The Maclaurin series representation of the functions G(r) and H(r) 
and their derivatives may be expressed according to the following 


scheme, 
G CG nn) abate) nea ee 
Hi H(r) + DH(r) + D°H(r) + D°H(r) + cee. 
DG DG(r) + D°G(r) + D°G(r) + D'G(r) + 2... 
on B) 3 hh WW V (3-28) 
DH DH(r) + D-H(r) + D-H(r) + D’H(r) +.... 
D°G D Zao + jf 36 (pr) + f Me cr) + f G(r) + .. 
5 = (3-29) 
DH D Bio + D-H(r) + D Te + f= Hr) + 4. 


where 


(3-30) 


Substitution of equations (3-3) through (3-30) into equation (2-9) 


converts the vorticity transport equation to the following form. 
3 I Vr 
r - D La(r) +d ae 1 r D-G(r) + D’G(r) + oe. 
M 3 ul ] 
le 0 pH ee : 3-11 D-H(r) a DeH(r) we 
/r* 


D Sale, + D°G(r) + pYG(r) 4. Saye i/e 
a a ad| 2 D aaa) “ D-H(r) . D!H(r) a 1/2 


BAe) re 4 D°G(r) + .e6 rn 
[0] ea eee Los) MGR & ooe 


1/r° 
1/r 
DG(r) + Boca + D°G(r) + pYc(r) oueeuews 
M Ver 
Ls) DH(r) + D-H(r) + DH(r) + DUA(r) +... 7 
4 


(r) + D°6(r) " DG(1) 

DG it ieG (re ct wae ir) =) ere ie 

Mi 

“tal fis ca od} * brtdiee 2247 


. 
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1/ri! 


1/r° 


36 (r) = pt 


OG 2 UG Geieamteliey) = 0 
moe nie 2 


G(r) + 1/2r 
+ DH(r) + D°H(r) + D-H(r) + DeH(r) +... 4) 1/6r 


1/2 
V/r° 
Ret ae (eC iaiea. aller 
* Po. lI, Pa Ce eet a(n) ee ea /2 
Chee) ee j r) + ré 
M 
ran : SF toa Th 
Whe 


Deke oes Meee. 
Yo] Done) + oe : “ae ] DH{r) + 20. 
Gia > eee eGo eames )) | yee 
¥[No-2| Mr) + DHE) + DoW) +... 1) 1/2 


Too ea 


Since equations (3-31) must be satisfied at all points in the flow, 


they must be satisfied in the limit as r approaches zero. A method of 
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determining the conditions to be satisfied at the singular point, r=0, 
is described below. 

Noting that as r approaches zero those terms in equations (3-31) 
containing r to the first power or greater may be neglected, these 
equations may be regrouped as a power series in r and expressed in the 


following abbreviated form, 


a ot 3 
c ts C 
L eel L ae 
. p°G(0) G(0} = 
D°H(0) H(0) 
D°G(0) DG(0) 
+ { (Cc + (Cc 
[os] p°H(0) L%} DE(O) 
\ D 


rs 


os 
O 
CO 
to 
ND 
Go 
vOvvxu—, 
© 
+ 
eS 
Lf 
ca 
vqOr 
GS 
i 
©) 
vx 
Wo 
WS) 
nO 
ee” 


where 


(Ch 8 ica Bitar (3-33) 


~HABP -AP(P - h) 


Re Re 
-p°p BP 
Re ke 


(co) = [4_3) * V5), (3-3)) 


nee) ees = 


re Re 
_pl 2B? 
re Re 
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fie) = (45.5) * 3) 7072). | (3-35) 


= | -2ABP -AP(P - )) 


Re ore 
_B°P 3BP 
eRe Re 


(o,] = Uo-2] -¥ Leal (3-36) 


= | 2ABQ AQ(P = 2) 
Reon: aeny Oo - | 
Re Re 
Cos} = Pt.) + Pe) + 072003) + O76) fMo a] (3-37) 


= -2AB(P + 3) ~13°(B" vel 
3Re 6Re 


-B°(P +3) 2B (P + 3) 
6Re She 


ce) = Pa) +P o-2) > ¥(f'427) , [%o-2) (3-38) 


= [oasc apg 
B°Q -2B0 
Cor} = (Mo) * Par) + O72) + 1/6) 3) + 728) [Hoy | 


= | AB (P + 8) -A(P + 8)(P + hh) 
OE 


Re euRe 
5 (3-39) 
-B° (P + 8) cK(P + 8) 
2Re OlRe 
[os] = oo) * Fred] + 072) [Mog] - 
¥(WN,_o) + Ly] + (1/2){No_o\) (3-0) 
= | ABQ NOK Be 2) 
, 
1 (Re : 30°) 5/30 3 s-) 
re ? Re 


(55) = Wo-o} - ¥B%-0 | (3-1) 
-lAB‘ Ac (* Q + | = AG 2)) 
Re 
i ia BO ; 
A (: {3 A \ 2B ) 2AB 


iN 


and 
P = Bo + 1 (3-2) 


Oe Bil ( -A \ + B& (3-13) 
Re 


As a consequence of equations (3-32), the following conditions must 


Beemet in the limit as r apvroaches zero. 


G(0) 0 
Lc,] a (3-hh) 
H(O) O 
, DG(0) {*} 
c 3-1 
[ 2 pees) (3-5) 


D°G(0) (0) 0 
Paco} * Mou het } a —— 
(D°c(0) DG(0) 0 
Cs T8.00) cH ; "| ial 
D'G(0) D°G(0) G(0) 0 
ic i + [ce] ” a + [c.) m {°\ (3-18) 
As discussed in Section II, the vorticity transport equation as ex- 
pressed by cauations (2-9) through (2-19) is, in general, a coupled set 
of differential equations. However, if Be0 these equations do uncouple. 
For this special case, the first of equations (2-9) becomes a homoze- 


neous fourth-order differential equation in H(r) and the second becomes 


a homogeneous second-order differential equation in G(r). Thus, the 


2h 


conditions expressed by equations (3-lh) through (3-8) may now be ex- 


amined for the special case B=O and the general case B>O. 


Inspection of equations (3-5) and (3-17) shows that they are 
identically satisfied for all values of the applicable odd-order 
derivatives of G(O) and H(O) for this case. Since the remaining 
equations uncouple for this case, the conditions on G(0) and H(0) may 
be studied incependently. 

Sue. ©) 

Sequential inspection of the first of equations (3-}l)), 
(3-6), and (3-8) verifies that there are three conditions to be en- 


forced. These three conditions ere as follows. 


H(O) = 0 (3-9) 

D°H(0) = 0 (3-50) 

D'y(0) = 0 (3-51) 
b. G(0) 


Although the second of equations (3-l) is identically 
satisfied for arbitrary G(0), sequential inspection of the second of 
equations (3-6) and (3-8) vields the following two conditions which 


must be enforced. 


G(O) = O (3-52) 
Bacon = 0 | (3-53) 
2. B>Oo 


With the exception of improbable special cases for which the 
determinants of any or all of the arrays [o,} through [Co} may be zero, 


the conditions to be met for this case are as follows. 
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G0) = H(0) = 0 (3-Sh) 


DG(O) = DH(O) = 0 (3-55) 
D°G(0) = D“H(O) = 0 (3-56) 
p°c(0) = D°H(0) = 0 (3-57) 
p''g(0) = D(0) = 0 | (3-58) 


The conditions expressed above represent a marked departure from 
the conditions vreviously thought to exist at r=O. Before the situation 
was properly understood, it was thought that the required conditions 
could be deduced from considerations of single-valuedness at this 
point. This approach is discussed in part two of Appendix G in Ref. 1 
Put must now be discarded as insufficient. In the interest of accuracy, 
the reader's attention is called to an error in equation (G-29) of that 
development. That equation should be corrected to read as follows. 


v, (0) ~ -iw, (0) (3-59) 
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IV. NUMERICAL METHODS 


Regrouping equations (2-9) allows expression of the equations to be 


solved in the following format. 


hy 


Tr a SP Ad GPs Bl 
M + | Pali 
bay di, Ls, ] ele zo, 


Dae _ (0c | the 
+ 0) + 
ne LS] = Yo) P 


Sane _ (0c (c 0 
= Oa 5a ey - + {X,| ' - ‘\ (h-1) 
The coefficient arrays are defined by equations (2-10) through (2-19). 

As discussed in Section III, the nature of this system of equations 
depends on the value of the complex PAE uan en wave number, B, For 
the special case B=0, the first of equations (lh-1) becomes a homoge- 
neous fourth-order differential equation in H(r) and the second becomes 
a homogeneous second-order differential equation in G(r). For the 
general case of B>O, equations (h-1) do not uncouple and must be sol- 
ved simultaneously. 

An additional influence of B is reflected in the character of the 
conditions at the singular point r=0. These conditions, as developed in 
Section III, are summarized below. 

For the special case B=0, the first of equations. (li-1) must simul-~ 


taneously satisfy the conditions 


H(O) = 0 (h-2) 
D°H(O) = 0 | (4-3) 
p!(0) = O (-h) 


and the second of equations (l-1) must similarly satisfy the following 


Cona: tions. 


of 


G(O) = 0 (4-5) 
Zz 
D°G(0) = 0 (4-6) 
For the general case where By0O, equations (l-1) must simultaneously 


satisfy the conditions below, 


G(O) = H(0) = 0 (-7) 
DG(O) = DH(O) = O (1-8) 
D°G(0) = D“H(O) = 0 (h-9) 
D°G(0) = D°H(O) = 0 (1-10) 
p'c(o) = DHH(O) = 0 | (1-11) 


The boundary conditions at the wall, r=1, are a consequence only of 
zero-velocity viscous effects at that point and thus do not vary with B. 


Maese cConaitions are as follows. 


G(1) = 0 (1-12) 
H(i) = 0 (1-13) 
DH(1) = 0 (-1h) 


Using the method of finite differences, the functions G(r) and H(r) 
may be approximated by a finite number of discrete, evenly spaced un- 
knowns. As shown in Figure li-1 below, the non-dimensionalized radius of 
the pipe may be divided into a one-dimensional computational mesh of 
uniform spacings consisting of n interior points, n+1 intervals, and n+2 
total points, including the boundary point at r=1 and the singular 
point at r=0O. The uniform spacing between each of the mesh points is 
§ as defined below. 


= 1/(n+1) (-15) 
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Figure 4-1 Finite Difference Mesh 
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By taking the Taylor series expansion of the function H(r) about 


the i 


th 


mesh point in terms of the values of this function at the mesh 


points 1+2, i+1i, i-1, and i-2 the second-order central difference ap- 


proximations of the derivatives of H(r) at the i 


have the following form. 


where 


eH. = 


ar 
we 
ae 
" 


Ss, 
ow 
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i+1 Hy _1)/28 


(Hay 


: a 
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Oop = 2a Teen eee 


Vile inlay oo Oo BL VE 


Ria lie 65 3505! 


My 


re 


th 


point are found to 


()-16) 
(hL-17) 
(4-18) 
(h-19) 


(1-20) 


The approximations for the derivatives of G(r) are determined in an 
identical manner. The error of equations (lj-16) through (4-19) is of 
the order of magnitude § °. The order of magnitude of this error may be 
reduced by expanding the method of derivation of these equations to in- 
clude more adjacent points on either side of the ae central = pount. 

It is important to note a discrepancy appearing in Ref. 1 concerning 
the development to this point. Careful comparison of Figure 3-1 in Ref. 
1 with Figure h-1 in this paper reveals that the direction of the 
labeling schemes used to depict the finite difference mesh has been re- 
versed. Note also that equations (3-2) in Ref. 1 compare exactly with 
equations (l-16) through (4-19) above. ‘When these equations are used in 
conjunction with the labeling scheme used in Figure 3-1 in Ref. 1, the 
signs of the odd-order derivative.approximations are reversed. This was 
amajor factor in producing the erroneous solutions obtained from the 
program presented in that paper. 

Substitution of the central difference approximations for the deri- 
vatives of the functions G and H in each of equations (l-1) results in 
a set of n linear, algebraic difference equations in terms of the un- 
known values for cach of these functions at each of the n interior 
points of the finite difference mesh depicted in Figure h-1. Since each 
of these equations is of the form of a linear combination of the ae 
central, mesh point and the two adjacent points on either side of this 
central point, this system of equations consists of banded coefficient 


arrays multiplying vectors containing the unknown values of the func- 


tions at each of the n interior points. By using this technique the 
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problem is converted to an eigenvalue problem of the general form 
Lxl{vt - yQl{v} = fot Ci) 
with the basic composition of the arrays [X] and (YJ and the vector {vt 


as illustrated in Figure h-2 below. 
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Figure h-2 Basic Composition of the Coefficient Arrays and 
the Vector of Unknowns. 


The exact composition of [X]} and [LY] depends upon the value of B. 
These arrays are established by the subroutine MSET2 in conjunction 
with the function subprograms CHM1£1 and CHM1E2, which compute the 
numerical value for each element in these arrays. The function sub- 
program CHM1E1 provides those coefficients required from the first of 
equations (li-1) and CHM1E2 provides those coefficients required from 
the second of these equations. For the special case B=0, all coeffi- 
cients contained in the upper right and lower left quadrants of Figure 


l-2 would be zero. 


ou 


Of particular interest are the difference equations whose central 
points are adjacent to the boundary at r=1, point n, and the singular 
point at r=O, point 1. To amplify this matter, consider the homoge- 
neous differential equation in H(r) resulting from setting B=O in 
equations (l-1). Use of equation. (l-19) to evaluate ‘the fourth deri- 


vative of this function at point n of the finite difference mesh re- 


veals that 
Hee 7s in 
Dew, = (H,, - WH, + 6H, - MHL, + H_p)/s'. (1-22) 
Similarly, for i=1 equation (l-19) becomes 
Din, - (Hy - LH, + 6H, - 4H, + H_,)/$". (1-23) 


Since Of r£1, equations (h-22) and (h-23) present apparent problems in 


that a is located at r=i+8 and H_, is located at r=-§, i.e. beyond 


2 


the allowable range of values for r. 


] 


These inconsistencies are resolved using the boundary conditions 
expressed by equations (lj-12) through (l-1h) and the conditions at r=0 
expressed by equations (h-2) through (lj-h). Using the labeling con- 
vention of Figure h-1 and equation (lj-13) it is easily confirmed that 

eC} = De Oey (j-2);) 


Comparison of equations (l-1h) and (l-16) yields the following relation- 


ships. 
DH(1) = Dees Ns - H)/2 = ¢ (-25) 
Snap 2 ly | (4-26) 


Equation (h-26) expresses the virtual point at r=1+§ in terms of the 
interior point n. Equation (h-22) may now be expressed exclusively in 
terms of boundary and interior points as follows (see depiction of 


equation at r=1-§ in Figure h-3 below). 
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Equation at r= 
Equation at r= 
Equation at r= 
Equation at r= 
Equation at r= 


Eouation at r=5§ 
Equation at r=l§ 
HauaLlon atl r=35 
Equation at r=28 
Equation at r=§ 


Central point o* the difference equation. 

Other points in the difference equation. 

Boundary or singular point where the value of the 
variable is zero. 

O Virtual point whase coefficient is combined with the 
coefficient of the point indicated by an arrow. 
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Figure h-3 Illustration of Basic Method Used to Include 
Virtual Points in Difference Equations. 
The virtual point at r=-§ can be expressed in terms of the interior 
point at r=§ in a similar manner. Expressed in terms of the labeling 


convention of Figure h-1, equation (l-2) odecomes 


Ba Os (h-27) 
Comparison of equations (})-3) and (4-17) shows that 
2 | an ; 
ie ot, = 2, + Ho Oe (4-28) 


Direct substitution of equation (-27) in equation (4-28) implies 
H me -H e ()-29) 
This leads to a more appropriate expression of equation (h-23) in terms 


of interior points as follows. 


ly ly 
D‘H, = (H, - WH, + 6H, - LH. - H,)/g (4-30) 
1 3 2 1 Ome 


Consider now the homogeneous second-order differential equation in 


G(r) for the special case B=0. Since the highest order derivative of 


Be. 


this function contained in equations (li-1) is second-order, inspection 
of equations (l-16) and (4-17) verifies the fact that no virtual points 
are required in the vicinities of the wall or the axis of symmetry. 

As previously mentioned, the general case of B>O requires solution 
of equations (h-1) as a system of coupled equations. Implementation of 
the conditions at r=0 for this problem as given by equations (h-7) 
through (li-11) requires special attention. 

The central difference approximation of the fourth derivative of H 
at r=§ requires a virtual point at r=-§. ixpression of this virtual 
point in terms of the interior point at r=$ is complicated by the re- 
quirement that the five conditions for the function H given by equations 
(4-7) through (lj-11) be satisfied simultaneously. Using equation (h-16) 
to satisfy equation (4-8) seems to imply that 

i ee (4-31) 
On the other hand, as previously shown in the developinent of equation 
(-29) 

H, = -H, . ()-32) 

This contradiction between equations (h-31) and (h-32) poses a pro- 
blem. However, the Maclaurin series approximation for the function H(r) 
and its Gerivatives in the vicinity of r=O is helpful in resolving this 


contradiction. These approximations are as shown below. 


2p acgoll a ee 6 ; 
H(r) = We Oe AoE = Aye —s + O(r’) (1-33) 
ree My 5 
BEM yoeaA ot Ar + Avo A ra + 7Al een Chesil 
a Ne 37 srl 
2 2 3 My 
Dence) = Aa + Alar tA, + Arete) (1-35) 
2 3 Ly O51 


3h 


D'H(r) = A, + Ayr + Agr’ + O(r°) (11-36) 
2 


De H(r) = Ay + Agr + OS ():-37) 


A comparison of equations (l-33) through ()-37) with equations ()-7) 
through (h-11) indicates that at r=0 
Ay =A, =A, = AZ = AL =O. —— (h-38) 


Substitution of equations (4-38) into equation (h-33) yields the fol- 


lowing sixth-order approximation for H(r) in the limit as r approaches 


ZETOs 

H(r) = here + o(r°) (4-39) 
Evaluation of equation (l-39) for r=+§ yields 

eco 

H(S) = Ag8 (11-100) 
share 

H(-$) = Ag(=8)° Gen) 

5! 


Inspection of equations (l-0) and (l-h1) quickly confirms that en- 
forcement of equations (l-7) through (l-11) requires that 

A) Shale (4-2) 
Note that the use of equation (h-39) as an approximation for the func- 
tion H(r) in equation (h-11) does not alter the second-order accuracy 
of the latter approximation. This may be conceptualized by a symbolic 


substitution of equation (4-39) in equation (1-11) as shown below. 


DYz, = (i, - LH, + 6, - Wy + ,) + 0(8°))/64 (4-113) 


(H, ~ LH, + 6H, - LH) + H_,)/s" * 0(82) (ebb) 
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With the formulation of the problem as discussed to this point 

accomplished, i.e. in the form of auaaeaen (+21), the method used 
for the remainder of the solution is summarized in the following steps: 
1.) Subroutine CDHTIN inverts the second coefficient array in equation 
(4-21), array [Y]. CDMTIN was obtained from the IPM Library routine 
CMTRIN by modifying it to make it applicable to double precision arrays. 
2;}:-Beth coefficient arrays, [X]} and [VY], are then premultiplied by 
‘ale Since multivlication of an array by its inverse invariably re- 
sults in the identity matrix, [I} , only the product ry fx) is com- 
puted using subroutine MULM. This converts the eigenvalue problem of 
equation (1-21) to the more conventional form 

(LZ3 - yf) (wy = fol. (4=-h5) 
where | 

U2) = fry") (4-6) 
3.) Since all the programs available for solving equations (l-h5) re- 
quire that the real and imaginary parts of the elements of [Z} be pre- 
sented in separate arrays, subroutine DSPLIT is called to accomplish 
it S 5 
4.) The eigenvalues of equations (h-l5) are computed by subroutines 
EHESSC and ELRH1C which are available through the International Mathe- 
Mavical and Statistical Library. These subroutines reduce the matrix 
[Z\ into the complex upper Hessenberg form and then solve for the 
eigenvalues. The results are then maaan back to the main program, R2. 
5.) The eigenvalues thus returned are listed on the computer output 


and plotted on the complex plane using subroutine PLOTP from the IBM 


Library. 
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V. RESULTS 


Owing to time constraints, the number of computer solutions obtained 
was sufficient only to confirm the general validity of the program pre- 
sented in this paper. For this reason, a detailed discussion of flow 
stability will not be undertaken here. However, the following brief 
@scussion should provide sufficient information for interpreting the 
data presented in this section. 

The characteristics of the flow which govern the stability are de- 
fined by the input variables A, B, and Reynolds ee A complete set 
of eigenvalues, ¥ , is obtained for each chosen combination of these 
variables. Recalling the form of the perturoation velocity vector po- 
tential, ‘I, as presented in Section II, 

Wi(x,r,8,t) = (e F(r) + é G(r) + e,H(r) )exp(Ax + BQ + Xt) Gear 
where 

mew + ¥, ; (5-2) 
it is readily seen that positive values for the real part, Sp. of the 
complex perturbation frequency represent an exponential growth rate in 
time and negative values of ¥, represent an exponential decay rate. ‘The 
algebraically greatest value of ¥, contained in the solution set is 
therefore used as the stability criterion. The corresponding eigenvalue 
is said to be the least stable root. The flow is Ween cs be 
stable, neutrally stable, or unstable with respect to stationary co- 
ordinates according to whether the least stable root is less than, equal 
to, or greater than zero, respectively. For a discussion of stability 


with respect to moving coordinates see part D of Section II in Ref. 1. 


of 


All computations were made using a half channel mesh containing 30 
interior points and with Ay held constant at a value of one. Tables I 


through V present the real part of the least stable root for each of 
the solutions obtained using the program developed in this ooo In 
general, these solutions are encouraging in that those errors which had 
previously caused solutions to exhibit decreasing stability with de- 
creasing Reynolds number seem to have been resolved. Results appear to 
be reasonable in that they exhibit stability characteristics which are 
qualitatively similar to those obtained for the plane flow problem as 
Seraizied in Section IV of Ref. 1. Specifically, for the range of Rey- 
nolds numbers investigated, all flows were stable for A.=0 while insta- 
bilities were found for Ap=-0.05. These results are significant in 


that previous investigations of fully developed pipe flow have never 


accounted for the instabilities that are known to exist. 
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